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Thermal or chaotic hght sources emit radiation characterized by a sUghtly enhanced probabihty of 
emitting photons in bunches, described by a zero-delay second-order correlation function g^"^^ (0) = 2. 
Here we explore photon-coincidence counting statistics of thermal cavities in the ultrastrong-coupling 
regime, where the atom-cavity coupling rate becomes comparable to the cavity resonance frequency. 
We find that, depending on the system temperature and coupling rate, thermal photons escaping 
the cavity can display very different statistical behaviors, characterised by second-order correlation 
functions approaching zero or greatly exceeding two. 

PACS numbers: 42.50.Pq, 42.50.Ar, 85.25.-j, 03.67.Lx 



Thermal radiation deserves a special place in modern 
physics. In the search for a solution to the discrepan- 
cies between the observed energy spectrum of thermal 
radiation and the predictions of classical theory, Planck 
was led to introduce the revolutionary concept of quanta 
[T]. Thermal or chaotic light sources emit radiation that 
is characterized by an enhanced probability of emitting 
photons in bunches [2 . In the course of the successful at- 
tempt at explaining such effect, Glauber established the 
basis of quantum optics O |4] . Even more recently the 
study of thermal emission has continued to provide amaz- 
ing results. A thermal light-emitting source is often pre- 
sented as a typical example of an incoherent source. How- 
ever, it has been shown that the field emitted by a ther- 
mal source made of a polar material is enhanced by more 
than four orders of magnitude and displays first-order 
coherence in the near- field zone [F. Moreover, by intro- 
ducing a periodic microstructure into such a polar mate- 
rial, a thermal infrared source can be fabricated that dis- 
plays first-order coherence over large distances [6 . While 
first-order coherence and spectral properties of thermal 
sources can be manipulated and tailored, their second- 
order coherence is known to be completely absent 0] re- 
sulting in the small bunching described by ^^^^(0) = 2. 

Here we investigate photon-coincidence counting 
statistics of thermal cavities in the ultrastrong-coupling 
regime, where the strength of the interaction (^r between 
an emitter and the cavity photons becomes comparable 
to the transition frequency of the emitter uj^ or the fre- 
quency of the cavity mode ujq. Cavities with ultrastrong 
light-matter interactions, are attracting interest both in 
semiconductor and superconducting systems, due to the 
possibility of manipulating the physical properties of the 
cavity quantum electrodynamic ground state [7W2l . 

The photon statistics of chaotic sources (like thermal 
cavities) and also of lasers can be explained classically. In 
contrast, strongly nonlinear photonic systems can emit 
photons one by one well separated in time from each 
other when excited coherently or operating very far from 
thermal equilibrium. For the systems considered so far. 



such scenario is however known to not persist when the 
coupled system is driven by thermal noise induced by 
reservoirs at a given temperature T. Indeed, the stan- 
dard quantum optics Master Equation (ME), generally 
used to study the dynamics of cavity QED systems [13] , 
predicts for such systems g^'^\^) = 2 independently of 
temperature and coupling strength. The interaction be- 
tween atoms and cavity photons is most often neglected 
when considering the coupling of this system with an 
environment. Recently it has however been shown that 
this simplification, which leads to the standard quantum 
optics ME, can give origin to unphysical effects in the ul- 
trastrong coupling regime [M]. Another key issue is the 
failure of standard quantum optical normal order corre- 
lation functions to describe photodetection experiments 
for such systems [151 US] • The theoretical treatment of 
this regime thus requires a description that goes beyond 
the standard techniques of quantum optics. 

By exploiting generalized correlation functions as in- 
troduced in [16] and a ME that fully takes into account 
the qubit-resonator coupling [14 , we investigate the 
photon-coincidence counting statistics of thermal sources 
for arbitrary light-matter coupling. The system that we 
investigate consists of a single mode resonator coupled 
to a two-level quantum emitter. Each subsystem is cou- 
pled to independent thermal baths of harmonic oscilla- 
tors at the same temperature T. We concentrate on the 
zero-detuning uoq = uo^ and low-temperature cases, where 
more striking deviations from standard results appear. 
We have corroborated the generality of our findings by 
confirming that similar results are found for one cavity 
mode coupled to multiple emitters or one emitter cou- 
pled to multiple cavity modes, see supplementary mate- 
rial [29]. 

A particularly well suited technology for such an ex- 
periment are superconducting circuits [lZl[T8| which have 
recently emerged as an excellent platform for microwave 
on-chip quantum-optics experiments and where second- 
order correlation function measurements for quantum 
[T9H2TI and low temperature thermal fields [22| have been 
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performed using quadrature amplitude detectors. 

Model - The cavity QED system that we explore con- 
sists of a single-mode resonator in interaction with a two 
level system (TLS). This interacting system can be de- 
scribed by the Rabi Hamiltonian (assuming = 1), 



H — ujoa^a 



+ ^(a + a'^)(cr" +cr+) (1) .g<2'(0) 



where uoq and uo^ are the cavity and TLS bare frequen- 
cies, g is the coupling strength, while a{a^) and <j~((T+) 
denote the annihilation (creation) operators for the cav- 
ity and TLS. Since we are interested in probing photon- 
coincidence statistics for arbitrary light-matter interac- 
tions where the contribution of the counter rotating 
terms cannot be neglected, we do not make use of the 
rotating wave approximation (RWA) [31 . Recently, it 
was shown that in the limit of very large coupling, the 
ground state and the first excited state of this system can 
become quasidegenerate (23H25] . 

In order to study thermal emission as well as the statis- 
tics of thermal photons we calculate the normal-order 
correlation functions of the output field. Standard nor- 
mal order correlation functions were recently shown to 
not correctly describe the emission properties and photon 
statistics of systems in the ultrastrong coupling regime 
[16]. For example, an application of the standard pro- 
cedure to the ultrastrong-coupling regime would predict 
an unphysical stream of output photons even for a zero- 
temperature system: {a^a)^^^ = Tr[a'''a p^^^] ^ 0. Fol- 
lowing [16 , we employ correlation functions for the out- 
put fields that are valid for an arbitrary coupling strength 



by expressing the cavity field X = Xo(a 



(Xo is 



the rms zero-point field-amplitude) in the atom-cavity 
dressed basis. In particular, the set-up we have in mind 
is equivalent to the emission of a thermalized black- 
box, whose output is coupled to the vacuum of a one- 
dimensional waveguide. In this case, by defining aout(^) 
and aj^'^it) as the output and input operators, the input- 
output relations become, a^^^{t) = aj^^{t) — zy^X+ 
[16 , where X+ is the positive frequency component of 
the cavity field X. The positive and negative frequency 
components {X~ = (X^)"*") can be derived by expanding 
the operators in the dressed state basis, namely the eigen- 
states \j) of H as in Eq. ([T]) ordered according to increas- 
ing eigenenergies Uj. Specifically, the time derivative of 
X+ can be expressed as X~^ = —'i^jkyj^kjXjk\j){k\ 
where Xj^ = {j\X\k) and Aj^j = Uk — ujj. According to 
these input-output relations and for input fields in vac- 
uum, the normalized second order correlation function 
for the output field can be expressed as 



(X-(t)X-(t + r)X+(t + r)X+(t)) 

t^oo (X-(t)X+(t))2 ' ^ ^ 



Results - The thermal-equilibrium zero-delay correla- 
tion function ^^^-'(0) can be directly calculated once the 
thermal equilibrium density-operator is known ^28^. For 
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FIG. 1: (color online) g^'^\0) calculated as function of the 
temperature and coupling strength. The results are obtained 
at resonance (cjq = <^x) and for a steady- state imposing ther- 
mal equilibrium between cavity and TLS, i.e. Ta = Tx. No- 
tably, at thermal equilibrium, statistical properties are inde- 
pendent of the damping rates. For comparison we show, the 
plane of g^'^\0) = 2 which indicates the value that would 
result from a conventional ME, where a RWA is performed. 



a system at thermal equilibrium, statistical properties are 
related to the density matrix of the canonical ensemble, 
that is the most general way to describe such thermalized 
interacting systems. Once H is in its diagonal form, the 
density matrix of the canonical ensemble reduces to. 



PT ^ ^ij, 



(3) 



where ej is the j-th eigenvalue of Sij the Kronecker 
delta function, and Z = exp(— e^/Zc^T) the parti- 
tion function. Figure 1 displays the thermal-equilibrium 
zero-delay correlation function ^^^-^(0) as a function of 
the effective coupling g/ooo and temperature calculated at 
zero detuning {ujq = uj^). The calculated ^'^^^(O) exhibits 
striking differences from the standard value of ^^^^ = 2. 
Of particular interest is the region with large effective 
coupling g/ujQ > 0.4 and low temperature ksT /uoq < 0.1 
where g^'^^ (0) becomes very small. Such highly nonclassi- 
cal behavior of thermal photons opens the perspective to- 
wards the realization of thermal sources of single photons 
in circuit QED. This anomalous behavior originates from 
the tendency of the interacting quantum system towards 
the vacuum degeneracy for large couplings. Specifically 
for increasing coupling the energy of the first excited state 
converges towards that of the ground state, while the 
other energy levels remain well separated from that dou- 
blet (see Fig. 2b). Hence, at sufficiently low temperature, 
only the first excited state is significantly populated by 
thermal noise {{uJi — uJo)/kT is non- negligible only for 
i = 1). At the onset of vacuum degeneracy, the ground 
and first excited state are quantum superpositions with 
multi-photon components [24 and at first sight one might 
expect to observe bunching effects. However such pho- 
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FIG. 2: (color online) a) Contourplot of ^*^^'*(0) calculated 
with the same parameters as Fig. [l] Here one can dis- 
tinguish four regions: i) milky, with 1.999 < ^^^^0) < 2 
corresponding to the standard thermal result, ii) gray, for 
1 < ^*^^'*(0) < 1.999, iii) blue, with the sub-Poissonian values 
^^^^(0) < 1 and iv) red, with a ^^^^0) > 2. The markers 
in these four regions are useful to characterise the different 
behaviors of g^'^\r) and their (g/uo, ksT/huJo) values are: 
diamond (0.1,0.2), triangle (0.2,0.1), square (0.5,0.07), circle 
(0.9, 0.15). b) Energy eingenvalues of H as in equation ([T]) as 
a function g/coo. 



tons are mostly virtual and the use of generalized nor- 
mal order correlation functions shows that the transition 
|1) |0) can only emit one physical photon at a time 
resulting into ^^^^ <C 1. 

For large effective couplings and higher temperatures 
g^'^\0) becomes larger than the standard value, ^''^^^(0) = 
2. Also this behaviour can be understood from the en- 
ergy levels in Fig. (2b). In particular at g/oJo — 0.45, a 
crossing between levels 2 and 3 (2 ^ 3) occurs. More- 
over the energy us lowers further for increasing g, so 
that, if temperature is not too low, thermal noise can 
populate the state |3), giving rise to an increase of the 
second order correlation function ^^^^(0) due to the cas- 
cade transitions |3) |1) |0). An overview of the 
behaviors of ^'^^^(0) is given in the contour plot of Fig. 
[2j where four different regions are displayed: i) a re- 
gion for small g/oJo (labelled by the diamond) where 
the standard thermal result ^^^^ ~ 2 is recovered (here 
1.999 < ^^^^(0) < 2) and which broadens for increasing 
temperatures; ii) a sub-Poissonian or nonclassical region 
with g^'^^ < 1 (square) in the ultrastrong coupling regime 
and for sufficiently low temperatures; iii) a region with an 
intermediate regime (triangle) with 1 < ^^^^(0) < 2; and 
iv) a super-Poissonian region (circle) beyond the stan- 
dard value g'^'^^ = 2 for very large coupling and higher 
temperatures. We note here that a calculation of ^''^^^(0) 
within the rotating wave approximation and using the 
standard quantum optics ME [13] would always yield the 
value ^^^-'(0) = 2, regardless of the coupling strength and 
the temperature. 



Although the steady state properties can be obtained 
from Eq. (|3|, a comprehensive description of the time 
dependent dynamics is possible only using the ME ap- 
proach. While zero-delay correlation function can be 
directly inferred from the thermal density operator, a 
description of the time dependent dynamics of the open 
quantum system is required to calculate the time-delayed 
second order correlation function g'^'^\r). A viable de- 
scription of system-bath interactions typically requires 
an expansion in the system bath coupling. A suitable 
way to perform this perturbative expansion consists in 
writing the Hamiltonian in the basis of its eigenstates 
In this way we obtain the following ME [Ml [26] 

p{t) = i[p{t),H] + CaP{t) + C,p{t), (4) 

where Ca and Cx are Liouvillian superoperators describ- 
ing the losses and the thermal feeding of the system. 
They read, C^pit) = E,,fe>, n(A,„ T,)P[|fc)(j|]p(t) + 
^j,k>j (l+n(A,„T,))D[|j)(fc|]p(t) for c = a.a' with 
V[0]p = \{20pO^ - pO^O - O^Op). The relaxation co- 
efficients = 27rdc{Akj)al{Akj)\Cjj^\'^ depend on the 
spectral density of the baths, dc{Ai^j)^ and the system- 
bath coupling strength, ac{Akj)^ at the respective transi- 
tion frequency Aj^j =00^ — 00 j as well as on the transition 
coefficients Cjk = —i{j\{c — c'^)\k) {c = a,a~). n{Akj,Tc) 
is the thermal population at frequency A^j and T the 
temperature. In all examples reported in this work, we 
consider a cavity that couples to the momentum quadra- 
tures of fields in one-dimensional output waveguides, as- 
suming that i) the spectral density dc{Akj) is constant 
and al{Akj) oc A^j ii) the cavity-TSL system is at the 
thermal equilibrium, i.e. = T^. Hence the relaxation 
coefficients reduce to T^^ = 7c (A^j/cJo) where 7c 

are the standard damping rates. These assumptions cor- 
respond to a cavity coupled to a one dimensional output 
waveguide as e.g. in circuit-QED. In the ME Q we ne- 
glect contributions of dephasing noise and Lamb shifts as 
they do not affect significantly our findings. A first re- 
sult worth mentioning is that the steady state solution of 
Eq. Q reproduces the thermal density matrix Eq. (|3|. 
This confirms that the present approach is able to de- 
scribe correctly the system thermalization, thus provid- 
ing a validation of Eq. ([3|. We then solve Eq. Q using 
as initial conditions the steady-state collapsed density 
matrix p{t) = X~^p{t)X~ and use p{t + r) to calculate 
g^'^\T) via the quantum regression theorem. 

Figure [3] displays g^'^\r) for temperatures and cou- 
plings corresponding to the reported markers in Fig. [2] 
and for 7a = 7cc = O.OIujq. It is worth to notice the 
appearance of oscillations in the triangle and diamond 
regions. Such oscillations arise from the superposition of 
the possible decay channels. In fact, for these values of 
the coupling strength, the separation between the energy 
levels is not high enough to suppress thermal occupation 
of higher excited states. Hence decays from upper energy 
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FIG. 3: (color online) g'''^\r) calculated for couplings and 
temperatures corresponding to the markers in Fig. |2] Here 
the damping rates are 7a = 7x = O.OIcjq. 



manifolds into lower excited states are possible. In this 
way, e.g. the third excited energy level can decay into the 
second or into the first, enabling a quantum superposi- 
tion of possible pathways and giving rise to this kind of 
oscillating behavior. This explanation is corroborated by 
the fact that the frequency of the oscillations corresponds 
exactly to A12, i.e. the energy difference between the first 
and the second excited state of H. Yet, in the square re- 
gion g^'^\T) is almost zero and it behaves as for a usual 
TLS, moreover the oscillations in r are suppressed by the 
low thermal noise and by the great separation in energy 
between the energy eigenvalues. Instead, the g^'^\r) in 
the dot region shows a slow dacay for large r which is 
in agreement with the extremely narrow linewidth of the 
transition between the first excited and the ground state 
(|1) |0)) (see Fig. [4|. In fact, owins; to the remark- 
able lifetime difference between the transitions |1) |0) 
and 1 3) the probability to measure a second pho- 

ton at time r is substantial as the first detected photon 
emerges from the decay of |3) into |1) and the second 
from the decay of |1) into |0). The latter is delayed due 
to the long lifetime of the transition |1) |0). For this 
reason, the g^'^\r) in the red region remains bunched for 
such a long time. Such behavior is a consequence of the 
spectral density of the bath that, being a linear function 
of the frequencies, tends to narrow the spectral linewidth 
of the lower resonances. 

Thermal emission is also characterized by its power 
spectrum, i.e. the Fourier transform of the two time cor- 
relation S{uj) (X \imt^oo2^J^{E-{t)E^{t^r))e'^^dr 
We exploit again the quantum regression theorem to cal- 
culate the relevant two-time correlation function. Fig. |4] 
shows calulations of S{uj) with the same parameters as 
used in Fig. [3] and the corresponding values of tempera- 
ture and couplings for the different markers as reported 
in the caption of Fig. [2] As expected the heights of 
the spectra increase for increasing temperature. More- 
over the resonances have i) different linewidths, as one 
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FIG. 4: (color online) Emission spectra calculated for cou- 
plings and temperatures corresponding to the markers in Fig. 
[2] and for 7^=7^= O.OIcjq- Spectra heights are normalized 
by the maximum value of the lowest temperature spectrum 
(square marker). 



can see from the definition of the damping rates and ii) 
different heights. The latter is mainly a consequence of 
the thermal feeding. For a fixed temperature, the ther- 
mal occupation of a spectral resonance (i.e. its height) is 
determined by n(cj, T) that increases as the frequency u 
of the resonance decreases. The presence of two different 
decay times in the example with the dot marker becomes 
apparent via the different linewidths of the resonances 
contributing to the signal, see Fig. |4^. 

Finally we observe that the effects described here could 
also be measured by populating the system with a qu- 
asithermal field distribution realized by mixing a fixed 
frequency microwave tone with noise sources of different 
bandwidths [21]. 
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FIG. 5: (color online) g^'^\0) calculated for two TLS coupled 
with a single cavity mode. 
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FIG. 6: (color online) g^'^\0) calculated for three TLS coupled 
to a single cavity mode. 

SUPPLEMENTARY MATERIAL 

In this section we report calculations of ^'^^^(0) for 
many TLS coupled to a single cavity mode, and for a 
single TLS coupled to two cavity modes. The latter 
case is of particular interest because of the unavoid- 
able coupling of the TLS to higher modes that is of- 
ten present in a real experimental setup. In each graph 
we plot ^*^^^(0) as function of the temperature and of 
coupling strength, imposing ujq = oj^ and = T^, 
i.e. at thermal equilibrium. We show moreover in all 
the graphs the plane ^^^^(0) = 2 to highlight the dif- 
ferences to the standard ME calculation which applies 
a RWA. The generalization of our model to many TLS 
is straightforward, and described by the Hamiltonian 
H = ujoa^a + E^x V+ar + (a + at) g^^) (ar + a+). 

Figs. |5) [6] and [7| show the feasibility to reach a robust 
quantum regime, characterized by the strong antibunch- 
ing, that persists even in presence of many TLSs. It is 
worth to notice the onset of a significant superbunching 
for higher coupling strenght. This is due to the degener- 
acy of the involved energy levels, that increase with the 
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number of TLSs. 
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FIG. 7: (color online) g^'^\0) calculated for four TLS coupled 
to a single cavity mode. 
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FIG. 8: (color online) g^^\0) calculated for a TLS coupled to 
two cavity modes. 



Fig. [sjin turn shows the thermal ^^^^(0) calculated for 
a single TLS coupled with two cavity modes as described 
by the Hamiltonian H = c^Q^^^aJai +cjQ^^a2a2 -\-(jOxa~^a~ + 
[g'^^^ai + a\) + g^'^\a2 + 4)](^" + ^^)- 1^ particular, 
for this latter calculation we used a second bare cavity 
mode with frequency cJq^^ = 2uJq'^ and a coupling with 
strength ^^^^ = 2g^^\ Such ratios of coupling strength 
appear for example in circuit QED experiments [8^. It 
is worth to notice the enhancement of the antibunching 
even though we consider a second bosonic mode. The 
explanation of such scenario is that the presence of the 
second cavity mode is compensated by the stronger cou- 
pling that effectively increase the nonlinarity of the whole 
system. 



